#include <math.h>
#include "fft_s.h"
#include "complex.h"

int factorization(int N, int *&f)
{
	int nmax, fmax, i, j, n;
	int *ftmp;
	nmax = (int)(log((double)N) / log(2.)); // fprintf(out,"%d ",nmax);fflush(out);
	ftmp = new int[nmax];
	n = 0;
	j = N;
	fmax = (int)sqrt((double)j);
	i = 2;
	while (i <= fmax)
	{ // fprintf(out,"%d %d %d\n",i,j,fmax);fflush(out);
		if (j % i == 0)
		{
			j /= i;
			fmax = (int)sqrt((double)j);
			ftmp[n] = i;
			n++;
		}
		else
		{
			if (i > 2)
				i += 2;
			else
				i++;
		}
	}
	ftmp[n] = j;
	n++;
	f = new int[n + 1];
	for (i = 0; i < n; i++)
		f[i] = ftmp[i];
	delete[] ftmp;
	return n;
} // factorization---------------------------------------------------------

fftdata::fftdata(int iN)
{
	double pi = 3.14159265358979323846;
	N = iN;
	n = factorization(N, f);
	int i, j, l, q, Nl;
	t1 = new complex[N];
	t2 = new complex[N];
	c = new complex **[n];
	cr = new complex **[n];
	for (j = 0; j < n; j++)
	{
		c[j] = new complex *[f[j]];
		cr[j] = new complex *[f[j]];
		for (l = 0; l < f[j]; l++)
		{
			c[j][l] = new complex[N];
			cr[j][l] = new complex[N];
		}
	}
	j = n - 1;
	Nl = N / f[j];
	for (i = 0; i < N; i++)
		for (l = 0; l < f[n - 1]; l++)
		{
			q = i / Nl;
			c[j][l][i] = exp_im(2 * pi * q * l / f[n - 1]) / N;
			cr[j][l][i] = exp_im(-2 * pi * q * l / f[n - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f[j];
		for (i = 0; i < N; i++)
		{
			q = i / Nl;
			for (l = 0; l < f[j]; l++)
			{
				c[j][l][i] = exp_im(2 * pi * q * l * Nl / N);
				cr[j][l][i] = exp_im(-2 * pi * q * l * Nl / N);
			}
		}
	}
} // fftdata

fftdata::~fftdata(void)
{
	int j, l;
	for (j = 0; j < n; j++)
	{
		for (l = 0; l < f[j]; l++)
		{
			delete[] c[j][l];
			delete[] cr[j][l];
		}
		delete[] c[j];
		delete[] cr[j];
	}
	delete[] c;
	delete[] cr;
	delete[] f;
	delete[] t1;
	delete[] t2;
} //~fftdata

fftdata2D::fftdata2D(int iN1, int iN2)
{
	double pi = 3.14159265358979323846;
	N1 = iN1;
	n1 = factorization(N1, f1);
	N2 = iN2;
	n2 = factorization(N2, f2);
	int i, j, l, q, Nl;
	t1 = new complex[N1 * N2];
	t2 = new complex[N1 * N2];
	c1 = new complex **[n1];
	cr1 = new complex **[n1];
	c2 = new complex **[n2];
	cr2 = new complex **[n2];
	for (j = 0; j < n1; j++)
	{
		c1[j] = new complex *[f1[j]];
		cr1[j] = new complex *[f1[j]];
		for (l = 0; l < f1[j]; l++)
		{
			c1[j][l] = new complex[N1];
			cr1[j][l] = new complex[N1];
		}
	}
	for (j = 0; j < n2; j++)
	{
		c2[j] = new complex *[f2[j]];
		cr2[j] = new complex *[f2[j]];
		for (l = 0; l < f2[j]; l++)
		{
			c2[j][l] = new complex[N2];
			cr2[j][l] = new complex[N2];
		}
	}

	j = n1 - 1;
	Nl = N1 / f1[j];
	for (i = 0; i < N1; i++)
		for (l = 0; l < f1[n1 - 1]; l++)
		{
			q = i / Nl;
			c1[j][l][i] = exp_im(2 * pi * q * l / f1[n1 - 1]) / N1;
			cr1[j][l][i] = exp_im(-2 * pi * q * l / f1[n1 - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f1[j];
		for (i = 0; i < N1; i++)
		{
			q = i / Nl;
			for (l = 0; l < f1[j]; l++)
			{
				c1[j][l][i] = exp_im(2 * pi * q * l * Nl / N1);
				cr1[j][l][i] = exp_im(-2 * pi * q * l * Nl / N1);
			}
		}
	}
	j = n2 - 1;
	Nl = N2 / f2[j];
	for (i = 0; i < N2; i++)
		for (l = 0; l < f2[n2 - 1]; l++)
		{
			q = i / Nl;
			c2[j][l][i] = exp_im(2 * pi * q * l / f2[n2 - 1]) / N2;
			cr2[j][l][i] = exp_im(-2 * pi * q * l / f2[n2 - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f2[j];
		for (i = 0; i < N2; i++)
		{
			q = i / Nl;
			for (l = 0; l < f2[j]; l++)
			{
				c2[j][l][i] = exp_im(2 * pi * q * l * Nl / N2);
				cr2[j][l][i] = exp_im(-2 * pi * q * l * Nl / N2);
			}
		}
	}

} // fftdata2D

fftdata3D::fftdata3D(int iN1, int iN2, int iN3)
{
	double pi = 3.14159265358979323846;
	N1 = iN1;
	n1 = factorization(N1, f1);
	N2 = iN2;
	n2 = factorization(N2, f2);
	N3 = iN3;
	n3 = factorization(N3, f3);
	int i, j, l, q, Nl;
	t1 = new complex[N1 * N2 * N3];
	t2 = new complex[N1 * N2 * N3];
	c1 = new complex **[n1];
	cr1 = new complex **[n1];
	c2 = new complex **[n2];
	cr2 = new complex **[n2];
	c3 = new complex **[n3];
	cr3 = new complex **[n3];
	for (j = 0; j < n1; j++)
	{
		c1[j] = new complex *[f1[j]];
		cr1[j] = new complex *[f1[j]];
		for (l = 0; l < f1[j]; l++)
		{
			c1[j][l] = new complex[N1];
			cr1[j][l] = new complex[N1];
		}
	}
	for (j = 0; j < n2; j++)
	{
		c2[j] = new complex *[f2[j]];
		cr2[j] = new complex *[f2[j]];
		for (l = 0; l < f2[j]; l++)
		{
			c2[j][l] = new complex[N2];
			cr2[j][l] = new complex[N2];
		}
	}
	for (j = 0; j < n3; j++)
	{
		c3[j] = new complex *[f3[j]];
		cr3[j] = new complex *[f3[j]];
		for (l = 0; l < f3[j]; l++)
		{
			c3[j][l] = new complex[N3];
			cr3[j][l] = new complex[N3];
		}
	}
	j = n1 - 1;
	Nl = N1 / f1[j];
	for (i = 0; i < N1; i++)
		for (l = 0; l < f1[n1 - 1]; l++)
		{
			q = i / Nl;
			c1[j][l][i] = exp_im(2 * pi * q * l / f1[n1 - 1]) / N1;
			cr1[j][l][i] = exp_im(-2 * pi * q * l / f1[n1 - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f1[j];
		for (i = 0; i < N1; i++)
		{
			q = i / Nl;
			for (l = 0; l < f1[j]; l++)
			{
				c1[j][l][i] = exp_im(2 * pi * q * l * Nl / N1);
				cr1[j][l][i] = exp_im(-2 * pi * q * l * Nl / N1);
			}
		}
	}
	j = n2 - 1;
	Nl = N2 / f2[j];
	for (i = 0; i < N2; i++)
		for (l = 0; l < f2[n2 - 1]; l++)
		{
			q = i / Nl;
			c2[j][l][i] = exp_im(2 * pi * q * l / f2[n2 - 1]) / N2;
			cr2[j][l][i] = exp_im(-2 * pi * q * l / f2[n2 - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f2[j];
		for (i = 0; i < N2; i++)
		{
			q = i / Nl;
			for (l = 0; l < f2[j]; l++)
			{
				c2[j][l][i] = exp_im(2 * pi * q * l * Nl / N2);
				cr2[j][l][i] = exp_im(-2 * pi * q * l * Nl / N2);
			}
		}
	}
	j = n3 - 1;
	Nl = N3 / f3[j];
	for (i = 0; i < N3; i++)
		for (l = 0; l < f3[n3 - 1]; l++)
		{
			q = i / Nl;
			c3[j][l][i] = exp_im(2 * pi * q * l / f3[n3 - 1]) / N3;
			cr3[j][l][i] = exp_im(-2 * pi * q * l / f3[n3 - 1]);
		}
	while (j > 0)
	{
		j--;
		Nl = Nl / f3[j];
		for (i = 0; i < N3; i++)
		{
			q = i / Nl;
			for (l = 0; l < f3[j]; l++)
			{
				c3[j][l][i] = exp_im(2 * pi * q * l * Nl / N3);
				cr3[j][l][i] = exp_im(-2 * pi * q * l * Nl / N3);
			}
		}
	}

} // fftdata3D

fftdata2D::~fftdata2D(void)
{
	int j, l;
	for (j = 0; j < n1; j++)
	{
		for (l = 0; l < f1[j]; l++)
		{
			delete[] c1[j][l];
			delete[] cr1[j][l];
		}
		delete[] c1[j];
		delete[] cr1[j];
	}
	for (j = 0; j < n2; j++)
	{
		for (l = 0; l < f2[j]; l++)
		{
			delete[] c2[j][l];
			delete[] cr2[j][l];
		}
		delete[] c2[j];
		delete[] cr2[j];
	}
	delete[] c1;
	delete[] cr1;
	delete[] f1;
	delete[] c2;
	delete[] cr2;
	delete[] f2;
	delete[] t1;
	delete[] t2;
} //~fftdata2D

fftdata3D::~fftdata3D(void)
{
	int j, l;
	for (j = 0; j < n1; j++)
	{
		for (l = 0; l < f1[j]; l++)
		{
			delete[] c1[j][l];
			delete[] cr1[j][l];
		}
		delete[] c1[j];
		delete[] cr1[j];
	}
	for (j = 0; j < n2; j++)
	{
		for (l = 0; l < f2[j]; l++)
		{
			delete[] c2[j][l];
			delete[] cr2[j][l];
		}
		delete[] c2[j];
		delete[] cr2[j];
	}
	for (j = 0; j < n3; j++)
	{
		for (l = 0; l < f3[j]; l++)
		{
			delete[] c3[j][l];
			delete[] cr3[j][l];
		}
		delete[] c3[j];
		delete[] cr3[j];
	}
	delete[] c1;
	delete[] cr1;
	delete[] f1;
	delete[] c2;
	delete[] cr2;
	delete[] f2;
	delete[] c3;
	delete[] cr3;
	delete[] f3;
	delete[] t1;
	delete[] t2;
} //~fftdata3D

void fft3D_rep(complex *in, complex *out, fftdata3D *data, int flag1, int flag2, int flag3)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1, m1, m2, n12;
	if (flag1 == 0 && flag2 == 0 && flag3 == 0)
	{
		for (int i = 0; i < (*data).N1 * (*data).N2 * (*data).N3; i++)
			out[i] = in[i];
	}

	if (flag1)
	{
		for (m1 = 0; m1 < (*data).N2; m1++)
			for (m2 = 0; m2 < (*data).N3; m2++)
			{
				temp1 = (*data).t1;
				temp2 = (*data).t2;
				j = (*data).n1 - 1;
				Nl = (*data).N1 / (*data).f1[j];
				Nl1 = 1;
				for (i = 0; i < (*data).N1; i++)
				{
					temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f1[j]; l++)
						temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
							in[k + l * Nl + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).c1[j][l][i];
				}
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0 && flag2 == 0 && flag3 == 0)
						temp1 = out;
					Nl = Nl / (*data).f1[j];
					Nl1 = Nl1 * (*data).f1[j + 1];
					for (i = 0; i < (*data).N1; i++)
					{
						temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f1[j]; l++)
							temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
								temp2[k + l * Nl + q1 * Nl * (*data).f1[j] + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).c1[j][l][i];
					}
				}
			}
	}

	if (flag2)
	{
		for (m1 = 0; m1 < (*data).N1; m1++)
			for (m2 = 0; m2 < (*data).N3; m2++)
			{
				if (flag1 && (*data).n1 % 2)
				{
					temp1 = (*data).t2;
					temp2 = (*data).t1;
				}
				else
				{
					temp1 = (*data).t1;
					if (flag1)
					{
						temp2 = (*data).t2;
					}
					else
					{
						temp2 = in;
					}
				}
				j = (*data).n2 - 1;
				Nl = (*data).N2 / (*data).f2[j];
				Nl1 = 1;
				for (i = 0; i < (*data).N2; i++)
				{
					temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f2[j]; l++)
						temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
							temp2[m1 + (k + l * Nl) * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).c2[j][l][i];
				}
				if (flag1 == 0)
					temp2 = (*data).t2;
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0 && flag3 == 0)
						temp1 = out;
					Nl = Nl / (*data).f2[j];
					Nl1 = Nl1 * (*data).f2[j + 1];
					for (i = 0; i < (*data).N2; i++)
					{
						temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f2[j]; l++)
							temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
								temp2[m1 + (k + l * Nl + q1 * Nl * (*data).f2[j]) * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).c2[j][l][i];
					}
				}
			}
	}

	if (flag3)
	{
		if (flag1)
			n12 = (*data).n1 % 2;
		else
			n12 = 0;
		if (flag2)
			n12 = (n12 + (*data).n2) % 2;
		for (m1 = 0; m1 < (*data).N1; m1++)
			for (m2 = 0; m2 < (*data).N2; m2++)
			{
				if ((flag1 || flag2) && n12)
				{
					temp1 = (*data).t2;
					temp2 = (*data).t1;
				}
				else
				{
					temp1 = (*data).t1;
					if (flag1 || flag2)
						temp2 = (*data).t2;
					else
						temp2 = in;
				}
				j = (*data).n3 - 1;
				Nl = (*data).N3 / (*data).f3[j];
				Nl1 = 1;
				if (j == 0)
					temp1 = out;
				for (i = 0; i < (*data).N3; i++)
				{
					temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f3[j]; l++)
						temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] +=
							temp2[m1 + m2 * (*data).N1 + (k + l * Nl) * (*data).N1 * (*data).N2] * (*data).c3[j][l][i];
				}
				if (flag1 == 0 && flag2 == 0)
					temp2 = (*data).t2;
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0)
						temp1 = out;
					Nl = Nl / (*data).f3[j];
					Nl1 = Nl1 * (*data).f3[j + 1];
					for (i = 0; i < (*data).N3; i++)
					{
						temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f3[j]; l++)
							temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] +=
								temp2[m1 + m2 * (*data).N1 + (k + l * Nl + q1 * Nl * (*data).f3[j]) * (*data).N1 * (*data).N2] * (*data).c3[j][l][i];
					}
				}
			}
	}
} // fft3D_rep----------------------------------------------------------

void fftr3D_rep(complex *in, complex *out, fftdata3D *data, int flag1, int flag2, int flag3)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1, m1, m2, n12;
	if (flag1 == 0 && flag2 == 0 && flag3 == 0)
	{
		for (int i = 0; i < (*data).N1 * (*data).N2 * (*data).N3; i++)
			out[i] = in[i];
	}

	if (flag1)
	{
		for (m1 = 0; m1 < (*data).N2; m1++)
			for (m2 = 0; m2 < (*data).N3; m2++)
			{
				temp1 = (*data).t1;
				temp2 = (*data).t2;
				j = (*data).n1 - 1;
				Nl = (*data).N1 / (*data).f1[j];
				Nl1 = 1;
				for (i = 0; i < (*data).N1; i++)
				{
					temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f1[j]; l++)
						temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
							in[k + l * Nl + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).cr1[j][l][i];
				}
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0 && flag2 == 0 && flag3 == 0)
						temp1 = out;
					Nl = Nl / (*data).f1[j];
					Nl1 = Nl1 * (*data).f1[j + 1];
					for (i = 0; i < (*data).N1; i++)
					{
						temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f1[j]; l++)
							temp1[i + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
								temp2[k + l * Nl + q1 * Nl * (*data).f1[j] + m1 * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).cr1[j][l][i];
					}
				}
			}
	}
	if (flag2)
	{
		for (m1 = 0; m1 < (*data).N1; m1++)
			for (m2 = 0; m2 < (*data).N3; m2++)
			{
				if (flag1 && (*data).n1 % 2)
				{
					temp1 = (*data).t2;
					temp2 = (*data).t1;
				}
				else
				{
					temp1 = (*data).t1;
					if (flag1)
					{
						temp2 = (*data).t2;
					}
					else
					{
						temp2 = in;
					}
				}
				j = (*data).n2 - 1;
				Nl = (*data).N2 / (*data).f2[j];
				Nl1 = 1;
				for (i = 0; i < (*data).N2; i++)
				{
					temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f2[j]; l++)
						temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
							temp2[m1 + (k + l * Nl) * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).cr2[j][l][i];
				}
				if (flag1 == 0)
					temp2 = (*data).t2;
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0 && flag3 == 0)
						temp1 = out;
					Nl = Nl / (*data).f2[j];
					Nl1 = Nl1 * (*data).f2[j + 1];
					for (i = 0; i < (*data).N2; i++)
					{
						temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f2[j]; l++)
							temp1[m1 + i * (*data).N1 + m2 * (*data).N1 * (*data).N2] +=
								temp2[m1 + (k + l * Nl + q1 * Nl * (*data).f2[j]) * (*data).N1 + m2 * (*data).N1 * (*data).N2] * (*data).cr2[j][l][i];
					}
				}
			}
	}

	if (flag3)
	{
		if (flag1)
			n12 = (*data).n1 % 2;
		else
			n12 = 0;
		if (flag2)
			n12 = (n12 + (*data).n2) % 2;
		for (m1 = 0; m1 < (*data).N1; m1++)
			for (m2 = 0; m2 < (*data).N2; m2++)
			{
				if ((flag1 || flag2) && n12)
				{
					temp1 = (*data).t2;
					temp2 = (*data).t1;
				}
				else
				{
					temp1 = (*data).t1;
					if (flag1 || flag2)
						temp2 = (*data).t2;
					else
						temp2 = in;
				}
				j = (*data).n3 - 1;
				Nl = (*data).N3 / (*data).f3[j];
				Nl1 = 1;
				if (j == 0)
					temp1 = out;
				for (i = 0; i < (*data).N3; i++)
				{
					temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] = 0;
					q = i / Nl;
					k = i % Nl;
					for (l = 0; l < (*data).f3[j]; l++)
						temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] +=
							temp2[m1 + m2 * (*data).N1 + (k + l * Nl) * (*data).N1 * (*data).N2] * (*data).cr3[j][l][i];
				}
				if (flag1 == 0 && flag2 == 0)
					temp2 = (*data).t2;
				while (j > 0)
				{
					j--;
					temp = temp1;
					temp1 = temp2;
					temp2 = temp;
					if (j == 0)
						temp1 = out;
					Nl = Nl / (*data).f3[j];
					Nl1 = Nl1 * (*data).f3[j + 1];
					for (i = 0; i < (*data).N3; i++)
					{
						temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] = 0;
						q = i / Nl;
						k = i % Nl;
						q1 = q % Nl1;
						for (l = 0; l < (*data).f3[j]; l++)
							temp1[m1 + m2 * (*data).N1 + i * (*data).N1 * (*data).N2] +=
								temp2[m1 + m2 * (*data).N1 + (k + l * Nl + q1 * Nl * (*data).f3[j]) * (*data).N1 * (*data).N2] * (*data).cr3[j][l][i];
					}
				}
			}
	}
} // fftr3D_rep----------------------------------------------------------

void fft2D_rep(complex *in, complex *out, fftdata2D *data, int flag1, int flag2)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1, m;
	if (flag1 == 0 && flag2 == 0)
	{
		for (int i = 0; i < (*data).N1 * (*data).N2; i++)
			out[i] = in[i];
	}

	if (flag1)
	{
		for (m = 0; m < (*data).N2; m++)
		{
			temp1 = (*data).t1;
			temp2 = (*data).t2;
			j = (*data).n1 - 1;
			Nl = (*data).N1 / (*data).f1[j];
			Nl1 = 1;
			for (i = 0; i < (*data).N1; i++)
			{
				temp1[i + m * (*data).N1] = 0;
				q = i / Nl;
				k = i % Nl;
				for (l = 0; l < (*data).f1[j]; l++)
					temp1[i + m * (*data).N1] +=
						in[k + l * Nl + m * (*data).N1] * (*data).c1[j][l][i];
			}
			while (j > 0)
			{
				j--;
				temp = temp1;
				temp1 = temp2;
				temp2 = temp;
				if (j == 0 && flag2 == 0)
					temp1 = out;
				Nl = Nl / (*data).f1[j];
				Nl1 = Nl1 * (*data).f1[j + 1];
				for (i = 0; i < (*data).N1; i++)
				{
					temp1[i + m * (*data).N1] = 0;
					q = i / Nl;
					k = i % Nl;
					q1 = q % Nl1;
					for (l = 0; l < (*data).f1[j]; l++)
						temp1[i + m * (*data).N1] +=
							temp2[k + l * Nl + q1 * Nl * (*data).f1[j] + m * (*data).N1] * (*data).c1[j][l][i];
				}
			}
		}
	}
	if (flag2)
	{
		for (m = 0; m < (*data).N1; m++)
		{
			if (flag1 && (*data).n1 % 2)
			{
				temp1 = (*data).t2;
				temp2 = (*data).t1;
			}
			else
			{
				temp1 = (*data).t1;
				if (flag1)
				{
					temp2 = (*data).t2;
				}
				else
				{
					temp2 = in;
				}
			}
			j = (*data).n2 - 1;
			Nl = (*data).N2 / (*data).f2[j];
			Nl1 = 1;
			if (j == 0)
				temp1 = out;
			for (i = 0; i < (*data).N2; i++)
			{
				temp1[m + i * (*data).N1] = 0;
				q = i / Nl;
				k = i % Nl;
				for (l = 0; l < (*data).f2[j]; l++)
					temp1[m + i * (*data).N1] +=
						temp2[m + (k + l * Nl) * (*data).N1] * (*data).c2[j][l][i];
			}
			while (j > 0)
			{
				j--;
				temp = temp1;
				temp1 = temp2;
				temp2 = temp;
				if (j == 0)
					temp1 = out;
				Nl = Nl / (*data).f2[j];
				Nl1 = Nl1 * (*data).f2[j + 1];
				for (i = 0; i < (*data).N2; i++)
				{
					temp1[m + i * (*data).N1] = 0;
					q = i / Nl;
					k = i % Nl;
					q1 = q % Nl1;
					for (l = 0; l < (*data).f2[j]; l++)
						temp1[m + i * (*data).N1] +=
							temp2[m + (k + l * Nl + q1 * Nl * (*data).f2[j]) * (*data).N1] * (*data).c2[j][l][i];
				}
			}
		}
	}

} // fft2D_rep----------------------------------------------------------

void fftr2D_rep(complex *in, complex *out, fftdata2D *data, int flag1, int flag2)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1, m;
	if (flag1 == 0 && flag2 == 0)
	{
		for (int i = 0; i < (*data).N1 * (*data).N2; i++)
			out[i] = in[i];
	}

	if (flag1)
	{
		for (m = 0; m < (*data).N2; m++)
		{
			temp1 = (*data).t1;
			temp2 = (*data).t2;
			j = (*data).n1 - 1;
			Nl = (*data).N1 / (*data).f1[j];
			Nl1 = 1;
			for (i = 0; i < (*data).N1; i++)
			{
				temp1[i + m * (*data).N1] = 0;
				q = i / Nl;
				k = i % Nl;
				for (l = 0; l < (*data).f1[j]; l++)
					temp1[i + m * (*data).N1] +=
						in[k + l * Nl + m * (*data).N1] * (*data).cr1[j][l][i];
			}
			while (j > 0)
			{
				j--;
				temp = temp1;
				temp1 = temp2;
				temp2 = temp;
				if (j == 0 && flag2 == 0)
					temp1 = out;
				Nl = Nl / (*data).f1[j];
				Nl1 = Nl1 * (*data).f1[j + 1];
				for (i = 0; i < (*data).N1; i++)
				{
					temp1[i + m * (*data).N1] = 0;
					q = i / Nl;
					k = i % Nl;
					q1 = q % Nl1;
					for (l = 0; l < (*data).f1[j]; l++)
						temp1[i + m * (*data).N1] +=
							temp2[k + l * Nl + q1 * Nl * (*data).f1[j] + m * (*data).N1] * (*data).cr1[j][l][i];
				}
			}
		}
	}
	if (flag2)
	{
		for (m = 0; m < (*data).N1; m++)
		{
			if (flag1 && (*data).n1 % 2)
			{
				temp1 = (*data).t2;
				temp2 = (*data).t1;
			}
			else
			{
				temp1 = (*data).t1;
				if (flag1)
				{
					temp2 = (*data).t2;
				}
				else
				{
					temp2 = in;
				}
			}
			j = (*data).n2 - 1;
			Nl = (*data).N2 / (*data).f2[j];
			Nl1 = 1;
			if (j == 0)
				temp1 = out;
			for (i = 0; i < (*data).N2; i++)
			{
				temp1[m + i * (*data).N1] = 0;
				q = i / Nl;
				k = i % Nl;
				for (l = 0; l < (*data).f2[j]; l++)
					temp1[m + i * (*data).N1] +=
						temp2[m + (k + l * Nl) * (*data).N1] * (*data).cr2[j][l][i];
			}
			while (j > 0)
			{
				j--;
				temp = temp1;
				temp1 = temp2;
				temp2 = temp;
				if (j == 0)
					temp1 = out;
				Nl = Nl / (*data).f2[j];
				Nl1 = Nl1 * (*data).f2[j + 1];
				for (i = 0; i < (*data).N2; i++)
				{
					temp1[m + i * (*data).N1] = 0;
					q = i / Nl;
					k = i % Nl;
					q1 = q % Nl1;
					for (l = 0; l < (*data).f2[j]; l++)
						temp1[m + i * (*data).N1] +=
							temp2[m + (k + l * Nl + q1 * Nl * (*data).f2[j]) * (*data).N1] * (*data).cr2[j][l][i];
				}
			}
		}
	}

} // fftr2D_rep----------------------------------------------------------

void fft_rep(complex *in, complex *out, fftdata *data)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1;

	if ((*data).n > 1)
		temp1 = (*data).t1;
	else
		temp1 = out;
	temp2 = (*data).t2;
	j = (*data).n - 1;
	Nl = (*data).N / (*data).f[j];
	Nl1 = 1;

	for (i = 0; i < (*data).N; i++)
	{
		temp1[i] = 0;
		q = i / Nl;
		k = i % Nl;
		for (l = 0; l < (*data).f[j]; l++)
			temp1[i] += in[k + l * Nl] * (*data).c[j][l][i];
	}

	while (j > 0)
	{
		j--;
		temp = temp1;
		temp1 = temp2;
		temp2 = temp;
		if (j == 0)
		{
			temp = temp1;
			temp1 = out;
		}
		Nl = Nl / (*data).f[j];
		Nl1 = Nl1 * (*data).f[j + 1];
		for (i = 0; i < (*data).N; i++)
		{
			temp1[i] = 0;
			q = i / Nl;
			k = i % Nl;
			q1 = q % Nl1;
			for (l = 0; l < (*data).f[j]; l++)
				temp1[i] += temp2[k + l * Nl + q1 * Nl * (*data).f[j]] * (*data).c[j][l][i];
		}
	}

} // fft_rep----------------------------------------------------------

void fftr_rep(complex *in, complex *out, fftdata *data)
{
	complex *temp1, *temp2, *temp;
	int i, j, k, l, q, q1, Nl, Nl1; // n=factorization(N,f);
	j = (*data).n - 1;
	Nl = (*data).N / (*data).f[j];
	Nl1 = 1;
	if ((*data).n > 1)
		temp1 = (*data).t1;
	else
		temp1 = out;
	temp2 = (*data).t2;
	for (i = 0; i < (*data).N; i++)
	{
		temp1[i] = 0;
		q = i / Nl;
		k = i % Nl;
		for (l = 0; l < (*data).f[j]; l++)
			temp1[i] += in[k + l * Nl] * (*data).cr[j][l][i];
	}
	while (j > 0)
	{
		j--;
		temp = temp1;
		temp1 = temp2;
		temp2 = temp;
		if (j == 0)
		{
			temp = temp1;
			temp1 = out;
		}
		Nl = Nl / (*data).f[j];
		Nl1 = Nl1 * (*data).f[j + 1];
		for (i = 0; i < (*data).N; i++)
		{
			temp1[i] = 0;
			q = i / Nl;
			k = i % Nl;
			q1 = q % Nl1;
			for (l = 0; l < (*data).f[j]; l++)
				temp1[i] += temp2[k + l * Nl + q1 * Nl * (*data).f[j]] * (*data).cr[j][l][i];
		}
	}

} // fftr_rep----------------------------------------------------------

void fft(complex *in, complex *out, int N)
{
	double pi = 3.14159265358979323846;
	complex *temp1, *temp2, *temp, *t1, *t2;
	int n, *f, i, j, k, l, q, q1, Nl, Nl1;
	n = factorization(N, f);
	if (n > 1)
		temp1 = new complex[N];
	else
		temp1 = out;
	if (n > 2)
		temp2 = new complex[N];
	t1 = temp1;
	t2 = temp2;
	j = n - 1;
	Nl = N / f[j];
	Nl1 = 1;
	for (i = 0; i < N; i++)
	{
		temp1[i] = 0;
		q = i / Nl;
		k = i % Nl;
		for (l = 0; l < f[j]; l++)
			temp1[i] += in[k + l * Nl] * exp_im(2 * pi * q * l / f[j]) / N;
	}
	while (j > 0)
	{
		j--;
		temp = temp1;
		temp1 = temp2;
		temp2 = temp;
		if (j == 0)
		{
			temp = temp1;
			temp1 = out;
		}
		Nl = Nl / f[j];
		Nl1 = Nl1 * f[j + 1];
		for (i = 0; i < N; i++)
		{
			temp1[i] = 0;
			q = i / Nl;
			k = i % Nl;
			q1 = q % Nl1;
			for (l = 0; l < f[j]; l++)
				temp1[i] += temp2[k + l * Nl + q1 * Nl * f[j]] * exp_im(2 * pi * q * l * Nl / N);
		}
	}

	if (n > 1)
		delete[] t1;
	if (n > 2)
		delete[] t2;
	delete[] f;
} // fft----------------------------------------------------------
void fftr(complex *in, complex *out, int N)
{
	double pi = 3.14159265358979323846;
	complex *temp1, *temp2, *temp, *t1, *t2;
	int n, *f, i, j, k, l, q, q1, Nl, Nl1;
	n = factorization(N, f);
	if (n > 1)
		temp1 = new complex[N];
	else
		temp1 = out;
	if (n > 2)
		temp2 = new complex[N];
	t1 = temp1;
	t2 = temp2;
	j = n - 1;
	Nl = N / f[j];
	Nl1 = 1;
	for (i = 0; i < N; i++)
	{
		temp1[i] = 0;
		q = i / Nl;
		k = i % Nl;
		for (l = 0; l < f[j]; l++)
			temp1[i] += in[k + l * Nl] * exp_im(-2 * pi * q * l / f[j]);
	}
	while (j > 0)
	{
		j--;
		temp = temp1;
		temp1 = temp2;
		temp2 = temp;
		if (j == 0)
		{
			temp = temp1;
			temp1 = out;
		}
		Nl = Nl / f[j];
		Nl1 = Nl1 * f[j + 1];
		for (i = 0; i < N; i++)
		{
			temp1[i] = 0;
			q = i / Nl;
			k = i % Nl;
			q1 = q % Nl1;
			for (l = 0; l < f[j]; l++)
				temp1[i] += temp2[k + l * Nl + q1 * Nl * f[j]] * exp_im(-2 * pi * q * l * Nl / N);
		}
	}

	if (n > 1)
		delete[] t1;
	if (n > 2)
		delete[] t2;
	delete[] f;
} // fftr----------------------------------------------------------
